home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / HQR.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  5KB  |  155 lines
  1. PROCEDURE hqr(VAR a: glnpnp; n: integer; VAR wr,wi: glnp);
  2. (* Programs using routine HQR must define the type
  3. TYPE
  4.    glnpnp = ARRAY [1..np,1..np] OF real;
  5.    glnp = ARRAY [1..np] OF real;
  6. where 'np by np' is the physical dimension of the matrix whose
  7. eigenvalues are to be found. *)
  8. LABEL 2,3,4;
  9. VAR
  10.    nn,m,l,k,j,its,i,mmin: integer;
  11.    z,y,x,w,v,u,t,s,r,q,p,anorm: real;
  12. FUNCTION sign(a,b: real): real;
  13.    BEGIN
  14.       IF (b < 0.0) THEN sign := -abs(a) ELSE sign := abs(a)
  15.    END;
  16. FUNCTION min(a,b: integer): integer;
  17.    BEGIN
  18.       IF (a < b) THEN min := a ELSE min := b
  19.    END;
  20. BEGIN
  21.    anorm := abs(a[1,1]);
  22.    FOR i := 2 TO n DO BEGIN
  23.       FOR j := i-1 TO n DO BEGIN
  24.          anorm := anorm+abs(a[i,j])
  25.       END
  26.    END;
  27.    nn := n;
  28.    t := 0.0;
  29.    WHILE (nn >= 1) DO BEGIN
  30.       its := 0;
  31. 2:      FOR l := nn DOWNTO 2 DO BEGIN
  32.          s := abs(a[l-1,l-1])+abs(a[l,l]);
  33.          IF (s = 0.0) THEN s := anorm;
  34.          IF ((abs(a[l,l-1])+s) = s) THEN GOTO 3
  35.       END;
  36.       l := 1;
  37. 3:      x := a[nn,nn];
  38.       IF (l = nn) THEN BEGIN
  39.          wr[nn] := x+t;
  40.          wi[nn] := 0.0;
  41.          nn := nn-1
  42.       END ELSE BEGIN
  43.          y := a[nn-1,nn-1];
  44.          w := a[nn,nn-1]*a[nn-1,nn];
  45.          IF (l = nn-1) THEN BEGIN
  46.             p := 0.5*(y-x);
  47.             q := sqr(p)+w;
  48.             z := sqrt(abs(q));
  49.             x := x+t;
  50.             IF (q >= 0.0) THEN BEGIN
  51.                z := p+sign(z,p);
  52.                wr[nn] := x+z;
  53.                wr[nn-1] := wr[nn];
  54.                IF (z <> 0.0) THEN wr[nn] := x-w/z;
  55.                wi[nn] := 0.0;
  56.                wi[nn-1] := 0.0
  57.             END ELSE BEGIN
  58.                wr[nn] := x+p;
  59.                wr[nn-1] := wr[nn];
  60.                wi[nn] := z;
  61.                wi[nn-1] := -z
  62.             END;
  63.             nn := nn-2
  64.          END ELSE BEGIN
  65.             IF (its = 30) THEN BEGIN
  66.                writeln('pause in routine HQR');
  67.                writeln('too many iterations'); readln
  68.             END;
  69.             IF ((its = 10) OR (its = 20)) THEN BEGIN
  70.                t := t+x;
  71.                FOR i := 1 TO nn DO BEGIN
  72.                   a[i,i] := a[i,i]-x
  73.                END;
  74.                s := abs(a[nn,nn-1])+abs(a[nn-1,nn-2]);
  75.                x := 0.75*s;
  76.                y := x;
  77.                w := -0.4375*sqr(s)
  78.             END;
  79.             its := its+1;
  80.             FOR m := nn-2 DOWNTO l DO BEGIN
  81.                z := a[m,m];
  82.                r := x-z;
  83.                s := y-z;
  84.                p := (r*s-w)/a[m+1,m]+a[m,m+1];
  85.                q := a[m+1,m+1]-z-r-s;
  86.                r := a[m+2,m+1];
  87.                s := abs(p)+abs(q)+abs(r);
  88.                p := p/s;
  89.                q := q/s;
  90.                r := r/s;
  91.                IF (m = l) THEN GOTO 4;
  92.                u := abs(a[m,m-1])*(abs(q)+abs(r));
  93.                v := abs(p)*(abs(a[m-1,m-1])+abs(z)
  94.                      +abs(a[m+1,m+1]));
  95.                IF ((u+v) = v) THEN GOTO 4
  96.             END;
  97. 4:            FOR i := m+2 TO nn DO BEGIN
  98.                a[i,i-2] := 0.0;
  99.                IF  (i <> (m+2)) THEN a[i,i-3] := 0.0
  100.             END;
  101.             FOR k := m TO nn-1 DO BEGIN
  102.                IF (k <> m) THEN BEGIN
  103.                   p := a[k,k-1];
  104.                   q := a[k+1,k-1];
  105.                   r := 0.0;
  106.                   IF (k <> (nn-1)) THEN
  107.                      r := a[k+2,k-1];
  108.                   x := abs(p)+abs(q)+abs(r);
  109.                   IF (x <> 0.0) THEN BEGIN
  110.                      p := p/x;
  111.                      q := q/x;
  112.                      r := r/x
  113.                   END
  114.                END;
  115.                s := sign(sqrt(sqr(p)+sqr(q)+sqr(r)),p);
  116.                IF (s <> 0.0) THEN BEGIN
  117.                   IF (k = m) THEN BEGIN
  118.                      IF (l <> m) THEN
  119.                        a[k,k-1] := -a[k,k-1];
  120.                   END ELSE BEGIN
  121.                      a[k,k-1] := -s*x
  122.                   END;
  123.                   p := p+s;
  124.                   x := p/s;
  125.                   y := q/s;
  126.                   z := r/s;
  127.                   q := q/p;
  128.                   r := r/p;
  129.                   FOR j := k TO nn DO BEGIN
  130.                      p := a[k,j]+q*a[k+1,j];
  131.                      IF (k <> (nn-1)) THEN BEGIN
  132.                         p := p+r*a[k+2,j];
  133.                         a[k+2,j] := a[k+2,j]-p*z
  134.                      END;
  135.                      a[k+1,j] := a[k+1,j]-p*y;
  136.                      a[k,j] := a[k,j]-p*x
  137.                   END;
  138.                   mmin := min(nn,k+3);
  139.                   FOR i := l TO mmin DO BEGIN
  140.                      p := x*a[i,k]+y*a[i,k+1];
  141.                      IF (k <> (nn-1)) THEN BEGIN
  142.                         p := p+z*a[i,k+2];
  143.                         a[i,k+2] := a[i,k+2]-p*r
  144.                      END;
  145.                      a[i,k+1] := a[i,k+1]-p*q;
  146.                      a[i,k] := a[i,k]-p
  147.                   END
  148.                END
  149.             END;
  150.             GOTO 2
  151.          END
  152.       END
  153.    END
  154. END;
  155.